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ABSTRACT 

We describe a three-dimensional, Godunov-type numerical magnetohydrodynamics (MHD) method 
designed for studying disk accretion to a rotating magnetized star in the general case where the star's 
rotation axis, its magnetic moment, and the normal to the disk all have different directions. The 
equations of ideal MHD are written in a reference frame rotating with the star with the z— axis aligned 
with the star's rotation axis. The numerical method uses a "cubed sphere" coordinate system which 
has advantages of Cartesian and spherical coordinate systems but does not have the singular axis of 
the spherical system. The grid is formed by a sequence of concentric spheres of radii Rj oc q J with 
j = 1..Nr and q = const > 1. The grid on the surface of the sphere consists of six sectors with the grid 
on each sector topologically equivalent to the equidistant grid on the face of a cube. The magnetic field 
is written as a dipole component plus deviations, and only the deviations are calculated. Simulation 
results are discussed for the funnel flows (FF) to a star with dipole moment /it at an angle O = 30° to 
the star's rotation axis l~i which is aligned with the normal to the disk. Results are given for different 
grids (Nr x N 2 ), from 26 x 15 2 (coarsest) to 50 x 29 2 (finest). We observe that the qualitative features 
of the accretion flows are very similar for the different grids, but the coarser grid is affected by numerical 
viscosity. We compare our 3D results for = with axisymmctric {2D), spherical coordinate system 
simulations of funnel flows (Romanova et al. 2002). Two important new 3D features are found in these 
simulations: (1) The funnel flow to the stellar surface is mainly in two streams which approach the star 
from opposite directions. (2) In the x — z cross section of the flow containing fi and 1~2, the funnel flow 
often takes the longer of the two possible paths along magnetic field lines to the surface of the star. 
A subsequent paper will give a detailed description of the method and results on 3D funnel flows at 
different inclination angles O. 
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1. INTRODUCTION 

There are many X-ray binary systems consisting of a 
rotating neutron star with dipole magnetic field which ac- 
cretes matter from the other star. Commonly, the accret- 
ing matter forms a disk as it approaches the neutron star. 
The neutron star's magnetic dipole axis fj, is not aligned 
with its rotation axis fl in the observed pulsating sources. 
Furthermore, in general /x and Jl are not aligned with the 
normal to the disk. That is, the dipole moment is inclined 
relative to the disk. In some systems, such as Her X-l this 
inclination is evident from the observations (e.g., Trumper 
et al. 1986). Models of magnetohydrodynamic (MHD) 
disk accretion to a rotating star with an aligned dipole 
magnetic field have been discussed by a number of authors 
(e.g., Ghosh & Lamb 1979; Camenzind 1990; Konigl 1991; 
Lovelace, Romanova, & Bisnovatyi-Kogan 1995, 1999; Os- 
triker & Shu 1995; Li & Wilson 1999; Koldoba et al. 2002). 
Analytical investigations of three dimensional (3D) accre- 
tion to an inclined magnetic rotator were done by Arons 
and Lea (1976a, b), Lipunov (1978a, b), Sharlemann (1978), 



Aly (1980), Lai (1999), and by Terquem and Papaloizou 
(2000). However, the essential aspects of this problem are 
three-dimensional, and these are not amenable to an ana- 
lytic approach. 

In this paper we describe a method developed specifi- 
cally for studying magnetohydrodynamic accretion to non- 
aligned magnetic rotators. The method is based on the 
"cubed sphere" coordinate system proposed by Sadourny 
(1972) and discussed by Ronchi, Iacono, & Paolucci (1996) 
(hereafter RIP96). We solve the equations of ideal mag- 
netohydrodynamics (MHD) in a reference frame rotating 
with the star using "cubed sphere" coordinates and us- 
ing a Godunov-type method. In contrast with RIP96, our 
simulations are in 31? rather than on the surface of the 
sphere, and furthermore we solve the equations of MHD. 

An important question is what grid resolutions are nec- 
essary for obtaining valid results for the problem of MHD 
disk accretion to a rotating star with magnetic moment 
mis-aligned with the star's rotation axis? Because 3D 
MHD simulations are very time-consuming, it is impor- 
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tant to find an optimal grid resolution. Here, we discuss 
results from simulations for the case where the star's mag- 
netic moment is inclined at an angle 9 = 30° to the star's 
rotation axis which is aligned with the normal to the disk. 
In §2 we describe our simulation methods. In §3 we sum- 
marize results of our 31? simulations of funnel flows to a ro- 
tating mis-aligned dipole, and we compare 3D simulation 
results for an aligned dipole with our earlier axisymmetric 
simulation results (Romanova et al. 2002). Conclusions of 
this work are given in §4. 

2. SIMULATION METHOD 

We consider a rotating magnetized star surrounded by 
an accretion disk and its corona. This problem is diffi- 
cult to treat numerically because the magnetic field varies 
strongly with distance from the star (~ 1/R?), and it is 
rapidly varying in the laboratory inertial reference frame. 
To minimize errors in calculating the magnetic force, the 
magnetic field B is decomposed into the "main" dipole 
component of the star, Bo, and the component, Bi, in- 
duced by currents in the disk and in the corona. Be- 
cause V x B = 0, the magnetic force is (J x B)/c = 
(V x Bi) x (Bo + Bi)/47r, which does not involve terms 
0(B 2 ) (Tanaka 1994; Powell ct al. 1999). 

Another difficulty in this problem is that the dipole mo- 
ment changes with time. It rotates with angular velocity fi 
so that the "main" field Bo also changes with time. Con- 
sequently in the induction equation there is a large term 
involving Bo- To overcome this difficulty we use a coor- 
dinate system rotating with angular velocity fi, in which 
the magnetic moment of the star \i and the "main" field 
B do not depend on time. 

2.1. MHD Equations in the Rotating Frame 

Let fi be angular velocity of rotation of the star. In the 
laboratory (inertial) reference frame the MHD equations 
are Dp/Dt + pV ■ u = 0, pDu/Dt = -Vp + (V x B) x 
B/4tt + pg, DS/Dt = , DB/Dt = (B • V)u - B(V • u) . 
Here, u is bulk velocity of the plasma, S the specific en- 
tropy, B the magnetic field, and D/Dt — d/dt + u • V the 
convective derivative in the inertial frame. 

The relations between the variables in the inertial frame 
and the reference frame rotating at the rate fi for a given 
fluid particle are u = v+fixR, B = B, where v is velocity 
of plasma in the rotating frame. The convective derivatives 
of scalar variables / do not change, Df/Dt = df /dt, where 
d/dt = d/dt + (v • V) is the convective derivative in the 
rotating frame. However, convective derivatives of vector 
variables transform as following: DF / Dt = dF/dt+flxF. 

Because V • (ft x R) = 0, the continuity equation is 
dp/ dt + p\I • v = . In the Euler equation in the rotating 
frame, there are two inertial terms added to the right hand 
side, the Coriolis force and the centrifugal force, pdv/dt = 
-Vp+(V x B) xB/47r+pg + 2 / 9 vxfi-p fix (fixR). In 
the induction equation (B • V)u = (B • V)v + (B • V)(fi x 
R) = (B ■ V)v + fi x B . Therefore in the rotating ref- 
erence frame, the induction equation has the same form, 
dH/dt = (B • V)v - B(V • v). 

Putting the equations in Eulerian form in the rotating 
frame, one finds that the new terms are the Coriolis and 
centrifugal forces in the Euler equation. The full set of 
equations is dp/dt + V • (pv) = 0, d(pv)/dt + V • T = 
pg + 2pvxfi-p fix (fix R), d(pS)/dt + V-(pSV) =0 , 



dH/dt = Vx(v x B), where T is the stress tensor with 
components T ik = p5 ik + pViV k + (B 2 6 ik /2 - BiB k ) / An. 
We do not include shocks in the present work as implied 
by the entropy conservation equation. 

2.2. The Grid 

The three dimensional grid consists of a set of concentric 
spheres of radii Rj with j = 1..Nr. The distribution of Rj 
is chosen to be inhomogeneous with Rj — R^q^ 1 where 
q =const. and R* is the radius of the numerical star. This 
choice implies ARj/Rj = q — 1. The grid on the surface 
of the sphere consists of six sectors with the grid on each 
sector topologically equivalent to the equidistant grid on 
the face of a cube. In each sector the grid of N x N cells 
is formed by the arcs of great circles separated by equal 
angles. For example, the +x sector is formed by the arcs 
of great circles going through the y and z— axes. 

This three-dimensional grid gives high spatial resolution 
close to the star which is important to our study of accre- 
tion to a rotating star with dipole magnetic field. Because 
ARj/Rj =const., the shape of the volume elements is ap- 
proximately cubical independent of R. 

The present study uses the scalings introduced by Ro- 
manova et al. (2002) where the inner radius of the disk and 
the beginning of the funnel flow (in the aligned, axisym- 
metric case) is at a radial distance i?o- Distances are mea- 
sured in units of Ro, velocities in units of «o = \JGM/ Rq, 
and time in units of Pq = 2ttRq/vq. With B the magnetic 
field at the surface of the star, we can define a reference 
density as po = B 2 /v 2 and a reference value for the mag- 
netic moment, ^o = B Tq. In this study the dimensionless 
radius of the numerical star is R m in/Ro = 0.35, and the 
outer radius of the simulation volume is R m nx/Ro ~ 4.8. 
We present results for two "cubed sphere" grids, one with 
N R x N 2 =50 x 29 2 cells in each of the six sectors (a total 
of « 2.5 x I0 5 cells), and the other with 26 x I5 2 cells 
per sector (a total of « 3.5 x 10 4 cells). For the first case 
q = 1.055 while for the lower resolution case q = 1.11. 

2.3. Boundary Conditions 

At the inner boundary R — i? m ; n - the "numerical star," 
- we assume "free" boundary conditions d/dR = for 
all variables. This boundary condition corresponds to ab- 
sorption of incoming matter so that there is no standoff 
shock. This boundary is treated as a rotating perfect con- 
ductor ft = flz. The flow velocity v was corrected such 
way that in the reference frame rotating with the star it 
is parallel to B at R — R m i n , that is, Bxv = at -R m i n . 
The boundary condition at i? m i n on the magnetic field has 
d(RB^)/dR = 0. At the outer boundary R = R max free 
boundary conditions are taken for all variables. 

2.4. Initial Conditions 

A star of mass M is located at the origin of the coordi- 
nate system. The initial magnetic field is the tilted dipole 
field of the star, B = 3(/x • R)R/i? 5 - fi/R 3 , where /x is 
tilted by an angle 9 from the z— axis (which is || fi). As 
initial conditions we set up a low-temperature Td, high- 
density pd disk with a high-temperature T c ^> Td, low- 
density p c <C Pd corona filling the remainder of the simu- 
lation region. The disk rotates with the same axis as the 
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Fig. 1. — Illustration of the "cubed sphere" grid which consists of six sectors. The +x sector is omitted in order to show the inner structure 
of the grid. The grid is inhomogencous in the R— direction in such a way that the cells are roughly cubical independent of R. 



Fig. 2. — Results of 3D MHD simulations of disk accretion to a rotating star with magnetic moment \JL tilted = 30° to the rotation axis 
f2. Two grids are shown, Nr x N 2 = 50 X 29 2 (the left-hand panels) and 26 X 15 2 (the right-hand panels) after t = 1.5Po, where Po is the 
period of rotation of the disk at Ro- The shading or color represents the density while the lines are magnetic field lines. The x — z cross section 
is the plane containing Si and fl. The x — y plane is the midplane of the disk at large distances from the star. Note that Nr is the number of 
grid cells in the radial direction while N X N is the number of cells on the surfaces of each of the six sectors of the "cubed sphere." 



Fig. 3. — Three dimensional view of the funnel flow to a rotating star with Si parallel to the z— axis and with the dipole moment \X tilted 
= 30° away from the z— axis in the x — z plane. The grid was Nr x TV 2 = 50 X 29 2 and the flow is shown at t = 2.5Po, where Po is the 
period of rotation of the disk at Ro- The thick arrow represents the star's magnetic moment. The red lines represent magnetic field lines. The 
yellow spiral line is a streamline. The nested blue and yellow lines are isodensity lines in the (x, z) plane. The central sphere represents the 
star which has a radius R t = 0.35i?o- The funnel flow close to the star is seen to be in two streams which approach the star from opposite 
directions. The density of the green surface shown is 0.35p0i where po is defined in §2.2. 



rotation of the star with angular velocity close to the Kep- 
lerian value to « Six ■ The disk extends inward to a radius 
Ro as discussed by Romanova et al. (2002) for the case of 
an aligned rotator. At this distance the ram pressure of 
the disk matter is of the order of the magnetic pressure of 
the dipole, p + pv 2 = B 2 /8tt. Initially, at any cylindrical 
radius r from the rotation axis, we rotate the corona and 
the disk at the same angular rate. This avoids a jump 
discontinuity of the angular velocity of the plasma at the 
boundary between the disk and the corona. Of course 
this distribution of u> in the corona leads to twisting of 
the dipole magnetic field lines. However, the twisting is 
smoothly distributed along the magnetic field lines, and it 
does not lead to fast evolution of the disk-corona system. 

2.5. Godunov-Type Finite-Difference Scheme 

All variables are evaluated at the centers of the cells. 
All vector variables are expressed in terms of their Carte- 
sian components. Finite difference equations are written 
for Cartesian components of vector variables. The finite 
difference scheme of Godunov's type has a form: 

l fl^l V+ £ s m T m = Q. (1) 

m— 1..6 

Here, U = {p, pv, B pS} is the "vector" of densities of 
conserved variables; T m is the "vector" of flux densities 
normal to the face "to" of the cube, s m is the area of the 
face "m" , V is the volume of the cell, Q is the intensity of 
sources in the cell, and At is the time step. 

To calculate the flux densities T m , an approximate Ric- 
mann solver is used, namely, an eight-wave Roe-type ap- 
proximate Riemann solver analogous to one described by 
Powell et al. (1999). The T's are represented as 

T m = X -{T' + T") - X -Y. ^\K\A a 1l a . (2) 

a 

Here, T' and T" are flux densities normal to the face 
"to" of the cell coming from neighboring cells, and v a are 
numerical coefficients which regulate the accuracy of the 
method. These are calculated using values of variables in 
the cells which are separated by this face. The A 's are ve- 
locities of the waves in the direction of the normal to the 
face; the A a 's are the wave amplitudes; and the 7£ a 's are 
the right-hand eigenvectors for the eight types of waves, 



a = (E, B, ±A, ±F, ±S) (see, e.g., Brio & Wu 1984; Powell 
et al. 1999). 

As a test of our 3D code, spherical Bondi accretion 
(Bondi 1952) was simulated for the case of a magnetic 
monopole field, B oc 1/R 2 with a grid 50 x ll 2 in each 
sector. The simulations shown good agreement with the 
radial dependence of the Bondi solution and with the pre- 
dicted accretion rate. 

3. SIMULATIONS OF FUNNEL FLOWS 

Here, we discuss simulation results for the case where the 
star's magnetic moment is inclined at an angle 9 = 30° to 
the rotation axis SI which is parallel to the normal of the 
disk. 

The magnetic moment of the star and the density of 
the disk are such that the ram pressure of the disk is ap- 
proximately equal to the magnetic pressure of the dipole 
at the inner radius of the disk Ro at t = 0. Thus, the 
funnel flow (FF) is expected to start from about this ra- 
dius, as observed in our axisymmetric simulations (Ro- 
manova et al. 2002). We rotate the star with angular 
velocity Q* = g{GM/Rl) 1 / 2 , with g = 0.19. The corota- 
tion radius of the star r cor is distance where the centrifu- 
gal force Sl 2 r equals the gravitational force GM/r 2 ; that is, 
r cor = (GM/ill) 1 / 3 . It follows that r cor = R /g 2/3 « 3i? - 

Some information about the 3D MHD flow can be ob- 
tained from the three orthogonal slices, an (x, z) slice at 
y = 0, an (y, z) slice at x = 0, and an (x,y) slice at 
z = 0. Figure 2 shows these cross-section for different 
grids. The most refined grid used in the present simula- 
tions has Nr x N 2 = 50 x 29 2 cells in each of the six sectors 
of the cubed sphere. The coarsest grid has 26 x 15 2 cells 
in each sector. Linking with a spherical coordinate sys- 
tem (R, 8, <fi), we have the correspondence for the number 
of cells, N R x N g x = 6N R x N 2 . We observed that 
the slices of the accretion flows are closely similar for the 
different grids. Namely, with decreasing distance, the disk 
ends and matter starts to flow out of the plane of the disk 
into the funnel flow at radii R w 2Ro in x— direction and 
at R w l.TRo in y— direction. In the (x,z) slice (which 
is the plane containing fi and SI) the funnel flow takes 
the longer of the two paths along the magnetic field. This 
we find to be a distinctive effect of the three-dimensional 
accretion flow to a mis-aligned dipole. Figure 3 shows a 
three dimensional view of the funnel flow. Close to the 
star the flow is in two streams which approach the star 
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from opposite sides. 

As a test of our code, we did 3D simulations of accre- 
tion to an aligned rotator (fi \\ ft) and analogous 2D ax- 
isymmetric simulations in spherical coordinates with the 
number of cells in the meridional plane Nr x Ng = 50 x 29. 
Results of the 2D and 3D simulations show qualitatively 
similar funnel flows at the same elapsed times. 

4. CONCLUSIONS 

We developed a new three dimensional MHD simulation 
code for studying the accretion flows to rotating stars with 
mis-aligned dipole magnetic fields. New three-dimensional 
features are found in the simulations of disk accretion to 
a rotating star with dipole moment fi inclined by an an- 



gle to the star's rotation axis. Specifically, in the x — z 
cross section of the flow containing fi and f2, the funnel 
flow takes the longer of the two possible paths along mag- 
netic field lines to the surface of the star. Furthermore, the 
funnel flow to the stellar surface is mainly in two streams 
which approach the star from opposite directions. A sub- 
sequent paper will give detailed analysis of the funnel flows 
at different inclination angles 6. 
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